package Math::Calculus::GaussElimination;

use strict;
use warnings;

use Moose;
use Moose::Util::TypeConstraints;
use Math::Cephes::Matrix qw/mat/;

our $VERSION = '0.01';

subtype 'ArrayNum' => as 'ArrayRef' => where {
    my $num = find_type_constraint('Num');
    ( $num->check($_) || return ) for @$_;
    return 1;
};

subtype 'ArrayArrayNum' => as 'ArrayRef' => where {
    my $array_num = find_type_constraint('ArrayNum');
    ( $array_num->check($_) || return ) for @$_;
    return 1;
};

has matrix => ( isa => 'ArrayArrayNum', is => 'rw', required => 1 );
has vector => ( isa => 'ArrayNum',      is => 'rw', required => 1 );

sub BUILD {
    my $self = shift;
    $self->validate;
}

sub validate {
    my $self = shift;
    my $mat  = mat( $self->matrix );
    my $a    = $mat->coef;
    for my $i ( 0 .. ( @$a - 2 ) ) {
	die "a.$i.$i cannot be 0" if $a->[$i]->[$i] == 0;
    }
    return 1;
}

sub iteration_estimate {
    my $self = shift;
    my $mat  = mat( $self->matrix );
    my $a    = $mat->coef;
    my $n    = @$a;
    return (4 * ( $n ** 3 ) + 9 * ( $n ** 2 ) - 7 * $n) / 6;
}

sub run {
    my $self = shift;
    my $mat = mat( $self->matrix );
    my ( $a, $b ) = ( $mat->coef, $self->vector );
    my $n = @$a;

    for ( my $k = 0 ; $k < $n - 1  ; $k++ ) {
	for ( my $i = $k + 1 ; $i < $n ; $i++ ) {
	    my $m = $a->[$i]->[$k] / $a->[$k]->[$k];
	    $a->[$i]->[$k] = 0;
	    for ( my $j = $k + 1 ; $j < $n ; $j++ ) {
		$a->[$i]->[$j] = $a->[$i]->[$j] - $m * $a->[$k][$j];
	    }
            $b->[$i] = $b->[$i] - $m * $b->[$k];
	}
    }

    my $x = $b;
    $x->[$n-1] = $b->[$n-1] / $a->[$n-1]->[$n-1];
    for (my $k = $n - 2 ; $k >= 0 ; $k-- ) {
        my $s = 0;
	for (my $j = $k + 1 ; $j < $n; $j++ ) {
	    $s += $a->[$k]->[$j] * $x->[$j];
	  }
	$x->[$k] = ( $b->[$k] - $s ) / $a->[$k]->[$k];
      }
    return $x;
}

1;
__END__


=head1 NAME

Math::Calculus::GaussElimination - Gauss Elimination Method

=head1 SYNOPSIS

  use Math::Calculus::GaussElimination;

=head1 DESCRIPTION

Perl extension for resolução numérica de sistemas lineares usando o gauss elimination method. Este método consiste en transformar o sistema linear original em um sistema linear equivalente com matriz dos coeficientes triangular superior, pois estes são de resolução imediata.

=head2 EXPORT

None by default.

=head1 SEE ALSO

Math::Cephes::Matrix

=head1 AUTHORS

Wallace Reis E<lt>wreis@cpan.orgE<gt>, Eden Cardim E<lt>edenc@cpan.orgE<gt>

=head1 COPYRIGHT AND LICENSE

Copyright (C) 2007 by Eden Cardim & Wallace Reis

This library is free software; you can redistribute it and/or modify
it under the same terms as Perl itself, either Perl version 5.8.8 or,
at your option, any later version of Perl 5 you may have available.

=cut
